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Abstract. - The probability P(a,N) that search algorithms for random Satisfiability prob- 
lems successfully find a solution is studied as a function of the ratio a of constraints per 
variable and the number N of variables. P is shown to be finite if a lies below an algorithm- 
dependent threshold «a, and exponentially small in N above. The critical behaviour is uni- 
versal for all algorithms based on the widely-used unitary propagation rule: P[(l + e) ot\, N] ~ 
exp( — iV 1 / 6 $(eA rl,/3 )). Exponents are related to the critical behaviour of random graphs, and 
the scaling function $ is exactly calculated through a mapping onto a diffusion-and-death 
problem. 



Introduction. The discovery of universality in phase transition phenomena was a major 
progress in modern condensed matter and statistical physics. The purpose of this letter is to 
point out that universality also takes place in computer science, more precisely, in compu- 
tational complexity theory. There, the goal is to understand whether a computational task 
consisting in processing a large number N of input data can be carried out in a time scaling 
only polynomially e.g. N 3 , and not exponentially e.g. 2 N , with N [1]. Depending on input 
data defining parameters, dynamical phase transitions between these two behaviours may take 
place [2-4]. We prove hereafter, for the case of the celebrated Satisfiability (SAT) problem [1], 
that the onset of complexity at criticality is universal in that it depends on some structural 
features of resolution algorithms and input data statistics only. 

Definitions of computational task and algorithm. In the random K-SAT problem [2], one 
wants to find a solution to a set of M — aN randomly drawn constraints (clauses) over a set 
of N Boolean variables Xi (i — 1 . . . N). Each constraint reads Zi x V Zi 2 V ... V Zi K) where 
V denotes the logical OR; Z£ is a variable X{ ( or its negation Xi e with equal probabilities 
(= |), and (ii,«2j • • • > ^k) is a i\-uplet of distinct integers unbiasedly drawn from the set of 
the (£) tf-uplets. We now study the K = 3 case, the smallest value for which the problem is 
NP-complete [1], and K > 4 later. 

Our algorithms start from tabula rasa, then iteratively assign variables to true (T) or 
false (F) according to two well-defined rules (specified below), and modify the constraints 
accordingly e.g. iiVi2Vi3 becomes X\ V23 if x^ — T [2,5]. At the end, no constraint is left (a 

© EDP Sciences 



2 



EUROPHYSICS LETTERS 



solution is found — success case) , or a contradiction is found (one variable previously assigned 
to, say, T is required to be F from modified constraints — failure case). The first assignment 
rule, UP (for unitary propagation) [5], is common to all algorithms: if a clause with a unique 
variable is produced at some stage of the procedure e.g. x\, then this variable is assigned 
to satisfy the clause e.g. xi = F. The second rule is a specific and arbitrary prescription 
taking over UP when it cannot be used i.e. in the absence of unique variable clause. In 
the simplest algorithm, referred to as R (random), the prescription consists in setting any 
unknown variable to T or F with prob. | independently of the remaining clauses [5]; more 
sophisticated prescriptions [6,7] will be studied later. 

Resolution procedures used in practical applications are based on the combination of the 
above algorithm and a backtracking principle [2]: in case of failure, the last variable assigned 
through the prescription (not through UP) is flipped, and the algorithm resumes from this 
stage. At the end, either a solution is found or all possible backtracks have failed, and a proof 
of the absence of solution is obtained. The resolution time typically scales as O(N) if a < a\ 
and exp O(N) if a > oa, where the threshold q;a depends on the algorithm. Intuition suggests 
and analyses prove [4,8] that this poly/exp crossover is due to the success/failure transition 
of the pure algorithm i.e. without backtracking. More precisely, oa can be identified with 
the ratio at which the probability P succ (a, N) of success of the pure algorithm vanishes as 
N — ► oo [5]. To understand the onset of complexity at qa, it is thus natural to analyze how 
-Psucc vanishes when the ratio a is kept close to its critical value and N increases. 

Analysis of the R algorithm. Each time a new variable is assigned some clauses are elim- 
inated, other are reduced or left unchanged. We thus characterize the set of clauses by its 
state (Ci(T),C2(T),C3(T)), where Cj is the number of j-clauses i.e. involving j variables 
(j = 1,2,3) and T is the number of assigned variables [4,5]. Consider a 3-clause left at 'time' 
T. When T — > T+l, the newly assigned variable has a probability 3/ (N—T) to appear (as is, or 
negated) in this 3-clause; if so the clause will be satisfied or reduced into a 2-clause (with equal 
prob. |). As a consequence the average change of C3 equals — 3Cs(T)/(N — T). In the large 
N limit, the density c%(t) = C%{tN)/N of 3-clauses becomes concentrated around its average 
value, solution of the ordinary differential equation dc^/dt — — 3c3/(l — t). A similar reasoning 
leads to dc2/dt — 3c3/2/(l — t) — 2c%/(l — t) for the density C2 of 2-clauses. Solving these equa- 
tions with initial conditions 03(0) = a, 02(0) = gives c${t) = a{\ — t) 3 ,C2{t) — §at(l — t) 2 
and the resolution trajectories of Fig.Q] [4]. 

The above evolution for C2, C3 is correct as long as no contradiction has emerged as a result 
of the production of two opposite 1-clauses e.g. x\ and x\. The probability Pn{C\]T) that 
the assignment of T variables has produced no contradiction and a set of constraints with C\ 
1-clauses obeys a Markovian evolution from T to T + 1 with a transition matrix [4] , 



H N [C{ «- C i; T,C 2 ] = 
+(l-Ic 1 )E- M Pr 1;si '° I c;-c 1+ i +Sl -, 2 J (1) 

si 

where Ic denotes the Kronecker function: Hp = f if C = 0, otherwise. Variables appearing 
in 0} are as follows: Sj (respectively rj) is the number of j-clauses which are satisfied (resp. 
reduced to j — 1 clauses) when the (T + l) th variable is assigned. These are stochastic 
variables drawn from multinomial distributions M.p ' x ' v = ( x C )p x+a (l — 2p) c ~ x ~~ y . Parameter 
Pj = j/2/(N — T) equals the probability that a j-clause contains the variable just assigned; 
7 1 ! = demands that no opposite 1-clauses and thus no contradiction are present. Equation 
p)l defines a random motion for a walker moving on the semi-infinite line C\ > with time- 
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Fig. 1 - Resolution trajectories for the R algorithm (bold lines, arrows indicate time direction, from 
a fraction t — of eliminated variables — right axis — to t = 1 — lower left corner where a solution 
is found). For initial ratio a < or, Ci stays bounded (success case). When a > art, Ci ~ N when 
the trajectory lies above the contradiction line 82 = 1, and the density c* is positive (failure case). 
At the critical ratio or(= |), the trajectory hits tangentially (black dot) the contradiction line. The 
critical region is defined by fluctuations ~ iV -1 / 3 for finite size N around these two lines (dotted and 
dashed lines respectively), and is crossed through assignment of a fraction At = N^ 1 ^ 6 of variables. 
Inset: C\ vs. T in the critical region for a particular run with N = 10 5 and a = qr. Reported 
scalings correspond to the largest components (S ~ iV 2 / 3 ). 

dependent and random (through C2) rates. The success/failure transition takes place when 
the average number 82 = c 2 /(l — <) of 1-clauses created from 2-clauses (right move) exceeds the 
number of 1-clauses (1 eliminated by UP each time a variable is assigned (left move) [5]. 

Successful regime (82 < l)- If elimination of 1-clauses is faster than creation, C\ stays 
bounded throughout the search process. On time scales 1 <C T <C N, C\ reaches equilibrium, 
with a distribution po(Ci, t) function of slow variables only i.e. C2, t. This, and the probability 
of success can be derived with the simple Ansatz Pn{Ci;T) = p (Ci,t) +p\(C\,t)/N + 
0(N~ 2 ) and sending — > 00. We find that po(Ci,t) ~ exp(— pC\) at large Ci, where 
p is the time-dependent positive root of p/ 82(f) = 1 — e~ p . As expected, p > and C\ 
is bounded as long as 8 2 < 1. At fixed ratio a, 8 2 reaches its maximum 8^ = |a along 
the resolution trajectory for a fraction tji = h of assigned variables; the transition takes 
thus place at ckr. = | [5]. The probability of success is given by P SUC c = Ylc Po(Ci,l) = 
exp[-^ — arctan(l/v / ?' — l)/2/\ / r — 1] with r = a^/a, and is shown in Fig.[2(top inset). Note 
that — lnP succ [aR(l — e)] ~ ? e -1 / 2 as a reaches its critical value an from below [9]. 

Failure regime (82 > !)■ When a > ckr, the resolution trajectory crosses the contradiction 
line 8 2 — 1 (Fig. . C\ then becomes of the order of A^ and each assignment has a finite 
probability to produce some contradiction; the probability of success is thus exponentially 
small with A^ [10]. Later the trajectory crosses the contradiction line back and C\ — 0(1) 
again. This behaviour is captured with the Ansatz Pjv(Ci;T) = exp(— Nu>(ci, t)) where u is 
the rate function associated to the large deviations of the density c\ of 1-clauses. ui fulfills a 



( 1 )More than one clause can be satisfied especially when C\ ~ N i.e. in the failure regime. 
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Fig. 2 - Scaling function $ (solid line) compared to numerical simulations for N = 1000 (empty), 
20000 (filled symbols) and algorithms R (A), GUC (□) and HL (CO- Error bars (with ~ 10 5 samples) 
are smaller than symbol size. Data for GUC and HL are rescaled horizontally and vertically, see text. 
Dotted lines serve as a guide for the eye. Bottom inset: — ln(P succ )/A rl//B vs. iV -1 / 6 at the critical 
thresholds aa, oguc — 3.003, «hl — 3.425. Linear fits (dotted lines) extrapolate to theoretical 
predictions for $(0) (available for R and GUC) on the left axis. Top inset: P SUC c vs. a showing the 
algorithm-dependent success/failure transition. Curves are analytical for R, GUC, and numerical for 
HL, KCNFS. 

first order partial differential equation [4,10], which can be explicitly solved for ratios slightly 
above the threshold i.e. a = aj\ (l + e). The coordinates c\,u>* of the minimum of lu i.e. the 
most likely trajectory are, at time t = ^(1 + Sy/e): c\ = 0, uj* = if s < — 1; c* = ~(1 — s 2 ) 2 e 2 , 

^ = i-i~i + i + TS^ 5/2 if N < !; c l = = A e5/2 if * > I- Thus, -lnP succ ~ ±e 5 ' 2 N 
to the lowest order in e. 

Critical regime (82 — 1). For N large but finite and a close to an, finite-size scaling [2] 
applies if 

- lnP succ ((1 + e) a R , N^j = N x $(e N e ) (2) 

for some regular function $. In the infinite size N limit, this expression should agree with the 
above results for the successful (e < 0) and failure (e > 0) regimes. Matching the powers in 
N and e, we find A — § = 0, $(2) - | l^l" 1 / 2 as x -> -00, and A + f = 1, - ^ x 5 / 2 as 
x — > +00. As a result, A = | and ^ = |. Figure (lower inset) shows the good agreement of 
numerical experiments performed at the critical point with the prediction — lnP succ ~ N 1 / 6 . 

Introduction of the oriented graph G representing 1- and 2-clauses allows us to understand 
the above scalings. G is made of 2 (N — T) vertices (one for each variable Xi and its negation 
Xi) , Ci marked vertices (one for each 1-clause z, ) , and 2 C2 edges (zi V Zj is represented by two 
oriented edges z\ — > Zj and Zj — * Zj) [11]. ^2 is simply the average (outgoing) degree of vertices 
in G. A step of UP corresponds to removing a marked vertex (and its attached outgoing 
edges), after having marked its descendents; steps are repeated until the connected component 
is entirely removed (no vertex is marked). Then a vertex is picked up according to the 
prescription and marked, and UP resumes. A contradiction arises when two 'opposite' vertices 
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i.e. associated to opposite variables are marked. The success/failure transition coincides with 
the percolation transition on G i.e. 8 2 = 1 as expected. From random graph theory [11,12], 
in the percolation critical window \82 — 1| ~ N~ x / 3 , the probability that a vertex belongs to 
a component of size S is Q(S) ~ S~ 3 / 2 , with a cut-off equal to the largest size, N 2 / 3 [13]. 
From Fig. ^ departure ratios a have to differ from cur by N 1 / 3 for resolution trajectories 
to fall into the critical window. Hence 8 — ^. The time spent by resolution trajectories 
in the critical window is At ~ yj\8 2 - 1| ~ A^ 1 / 6 , corresponding to AT = N At ~ TV 5 / 6 
eliminated variables. Let Si, S2, ■ ■ ■ , Sj be the sizes of components eliminated by UP in the 
critical window; we have J - AT/ JdSQ{S)S ~ N 1 / 2 . During the j th elimination, the 
number of marked vertices 'freely' diffuses, and reaches C\ ~ \[S] (Inset of Fig. [IJ. The 

probability that no contradiction occurs is [(1 — q) Cl } Sj ~ e~ s i' ^ N where q ~ jj is the 
probability that a marked vertex is 'opposite' to the one eliminated by UP. Thus — lnP succ ~ 
J J dS Q(S) S 3 / 2 /N ~ N 1 / 6 , giving A = ~. Notice that, while the average component size is 
S ~ A^ 1 / 3 (and thus P^(C\ — 0) ~ Af -1 / 3 ), the value of A is due to the largest components 
with S ~ N 2 / 3 i.e. G\ ~ A^ 1 / 3 marked vertices. The distribution of c = C1/N 1 / 3 is calculated 
below. 

Scaling function. To calculate $ in (J21, we magnify the critical region in Fig. ^ an d 
consider ratios a — |(1 + eo N~ e ) and times t = | (1 + to N~ T ), where 6> and r are scaling 
exponents to be determined. We then decompose the probability of having C\ clauses in the 
critical region into the product Pjv(Ci, T) — exp[— N x ip(to)] x F(Ci,to); the first term is the 
probability that no contradiction has been found up to 'time' to, and F is the (normalized) 
probability distribution of 1-clauses. Clearly <p(to — > —00) = since the probability that 
the search process has ended is not vanishingly small before the trajectory enters the critical 
region (Fig.^J. We make the Ansatz F(C\) = N^ 1 f(c = CiN~ 7 ) where / is the probability 
distribution of the rescaled number c of unit-clauses (Fig.[3Jl. Last of all, the probability that 
a variable is set through a free choice and not UP, P^{C\ =0), is assumed to scale as fo/N 10 . 

The evolution equation for Pn based on matrix Q imposes S = 7 = 7o = ^,A = r= | 
in agreement with the above scaling arguments ( 2 ). In addition, we find dip/dto — c(to), the 
average value of c with distribution / solution of 

1 d 2 f df 

1 «o^ + (c- C )/ = (3) 



2 5c 2 u d 



c 







with vq = t$ — €q. Boundary conditions are (d c f + Wo/)|c=o — (reflecting wall) and /( 
/(0)/2. The diffusion term in J2J reflects the Gaussian stochastic nature of 2- to 1-clauses 
reductions, the drift term favors small (respectively large) values of the density c when vo > 
(resp. wo < 0) — corresponding to 82 < 1 and 82 > 1 respectively — and the third term 
expresses the relative death-rate of search processes with respect to the average rate c (the 
higher c, the more likely it is to encounter a contradiction)( 3 ). The solution of differential 
equation © reads /(c) oc exp(— v c) Ai[y / 2c+ z(vo)], where Ai is the Airy function, z(x) is 
the inverse function of x(z) = — \[2 Ai'(z)/Ai(z). Distribution / is shown on Fig.|3for several 
values of the drift. Positive (respectively negative) Vq correspond to trajectories below (resp. 
above) the contradiction line (Fig-QJ, with distributions / peaked around c = (resp. c > 0). 

The probability of success remains unchanged (to the leading order in N) once the tra- 
jectory exits from the critical region, thus — lnP succ /A^ A = <p(to — * +00). This proves the 

( 2 ) These results are not affected by the presence of Gaussian fluctuations ~ N^ 1 / 2 in C2(i). 

( 3 ) Notice that lacks any time derivative since / reaches stationarity on a much shorter time-scale (TV -1 / 3 ) 
than the one relevant for P S ucc 
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number C, of 1 -clauses 



Fig. 3 - Histograms of the numbers Ci of 1-clauses for sizes N = 10 2 (□), 10 3 (O) and 10 4 (A) 
right at criticality (time t = | and ratio a = §). Note the sharp increase of the probability around 
Ci = 0, in quantitative agreement with the theoretical prediction fo = /(0)/2. Insets: theoretical 
distributions / for c = Ci/N 1 ^ 3 at criticality (A, same data as main figure), and for drifts vo = 0.5 
(B, dotted), vq = —0.7 (B, dashed curve) compared to numerics. 

existence of the scaling function defined in J2J , with 

^ = iCtB=J x2 -^^ (4) 

Scaling function $ is plotted and compared to numerics in Fig. [3 It is an easy check that all 
the results related to the successful and failure regimes e.g. the values c*,u* listed above for 
finite e and N — > oo are found back when eo — > ±°o respectively. 

Universality. The critical point of R, or any algorithm A that implements UP e.g. GUC [5] 
or HL [6] where variables are chosen to satisfy 2-clauses or according to their occurrences 
respectively, is reached when the resolution trajectory is, for some time t\, tangent to the 
82 = 1 line: ^(^A + Ai) — 1 = b(At) n with n > 2 even integer and b determined from 
the derivatives of the density C2 of 2-clauses at £a- The tangency condition reflects that the 
creation of new edges in G, from the reduction of 3- into 2-clauses, precisely compensates the 
elimination of edges by UP. Although the resolution trajectory strongly depends on A, the 
critical behaviour depends on n, b only, and is thus universal. 

The value of n is 2 for R, GUC, HL and generic algorithms A. Therefore, 6 = = | 
independently of A; the scaling function is <J>a(£o) = r % ®{ r A £o) where the ta's are functions 
of b e.g. rg uc ~ 0.9902, r^ uc ~ 1.7182. Fig. illustrates that GUC and HL fall in this 
UP universality class. Numerical investigations suggest that more complex algorithms as 
KCNFS [7] do, too. This universality class is robust against any change, either induced by 
algorithms or present in the input data distribution, in the degree sequence of the clauses graph 
G, a consequence of the robustness of the critical component size distribution Q(S) [13]. 

Higher values of n(> 4) are exceptionally found for finely-tuned input data statistics e.g. 
with clauses of different lengths £(< K) and appropriate ratios cti (so far we have restricted 
to K — 3 with «3 = a). If the ratios at the critical point 62 — 1 are such that the reduction 
of (£ + 1)- to ^-clauses compensates the disappearance of ^-clauses for all 2 < £ < K, then the 
resolution trajectory will stay longer in the critical region, making P succ decrease. The precise 
condition is a e = 2 e ~ 1 /l/(£ - 1) at criticality [14], leading to A = with n = K —I ( 2 ). 

It would be interesting to extend our study to structured input data distributions e.g. 
leading to clause graphs G embedded in finite-dimensional spaces, and possibly to 9 ^ ^. In 
this context, developing renormalization tools to capture the critical behaviour of algorithms 
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would be the natural yet apparently difficult next step. It would also be worth to study uni- 
versality for other types of algorithms e.g. local search procedures [15], or other computational 
tasks [3] e.g. graph coloring [9], where poly/exp transitions take place. 
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